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This paper is concerned exclusively with axisymmetric spacetimes. We want to develop reductions 
of Einstein's equations which are suitable for numerical evolutions. We first make a Kaluza-Klein 
type dimensional reduction followed by an ADM reduction on the Lorentzian 3-space, the (2+1) + 1 
formalism. We include also the Z4 extension of Einstein's equations adapted to this formalism. Our 
gauge choice is based on a generalized harmonic gauge condition. We consider vacuum and perfect 
fluid sources. 

We use these ingredients to construct a strongly hyperbolic first-order evolution system and 
exhibit its characteristic structure. This enables us to construct constraint-preserving stable outer 
boundary conditions. We use cylindrical polar coordinates and so we provide a careful discussion of 
the coordinate singularity on axis. By choosing our dependent variables appropriately we are able 
to produce an evolution system in which each and every term is manifestly regular on axis. 
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I. INTRODUCTION 
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In the early days of numerical relativity attention focussed on spherically symmetric spacetimes. The next simplest 
subject would have been axisymmetric spacetimes, but very little progress was made in this area because, in polar 
coordinates adapted to the symmetry, there is a coordinate singularity on axis. Many attempts to deal with this proved 
unsuccessful, most evolutions became unstable and attention quickly turned, with the rapidly increasing capacity and 
speed of hardware, to spacetimes without symmetries. For a comprehensive review see e.g., |l3j . 

However axisymmetric situations arise frequently in astrophysics and deserve further study. The first question is, 
5h ■ obviously, why were the evolutions unstable? Suppose we assume that in a neighbourhood of the axis the spacetime 
jpP" is regular so that we can introduce locally cartesian coordinates (t,z,x,y) with the axis given by x = y = 0. Then 
J> , the Killing vector is £ = —yd/dx + xd/dy, and we say that M is axisymmetric if C^M — 0. We say that a tensor 
field M is regular on axis if its components with respect to this chart can be developed as a Taylor series in x and 
y in a neighbourhood of the axis 1 . We may introduce polar coordinates r = y/ x 2 + y 2 and tp — arctan(y/a;), so that 
the Killing vector is £ = (d/dtp). Because the transformation between cartesian coordinates and polar coordinates 
(t,z,r,cp) is singular on axis 2 , axisymmetric tensor fields which are regular on axis may take strange forms when 
expressed in polar coordinates. For example if M a p is symmetric and has these properties then 



M af3 



( A B rD r 2 F \ 

B C rE r 2 G 

rD rE H + r 2 J r 3 K 

\ r 2 F r 2 G r 3 K r 2 (H - r 2 J) J 



(1) 



where A, B, . . . , K are functions of t, z and r 2 . (See appendix for this and related results.) In particular both 
the metric and its first time-derivative must take this form. If the evolution scheme fails to preserve precisely the 
indicated r-dependencies then M becomes irregular on axis and instability is inevitable. 



'Electronic address: o.rinne@damtp.cam.ac. 
^Electronic address: j.m. stewart@damtp.cam.ac.uk 

1 For numerical purposes analytic and regular are synonymous. 

2 Note that the cartesian expression for the Killing vector £ vanishes on axis, whereas the polar form does not. 
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Clearly the problems lie mainly in the last row and column of (JIJ. In section[H]we outline a reduction of Einstein's 
equations in spacetime M. to a 3-dimensional Lorentzian manifold M with chart (i, z, r), effectively removing the last 
row/column. (The relation between M rr and M„„ still needs careful treatment.) Given a solution of the reduced 
equations in Af one can recover the solution in AA. However there is no need to do this. All axisymmetric quantities 
of physical interest in M have their counterparts in Af, and so one might as well work in Af. Because that manifold 
is Lorentzian we can perform an ADM decomposition. The details of this (2 + 1) + 1 scheme, with an arbitrary 
axisymmetric matter source, are given in section ITT1 

In the (2 + 1) + 1 system we have a spatial foliation with lapse function a, shift vector i] A , and spatial 2-metric 
Hab- (Here A, B, . . . range over r and z.) Since every 2-metric is conformally flat it is tempting to set H rz = 0, 
H rr = H zz = if 4 say. The constraint equations imply elliptic equations for ip and the shift vector components, 
and a quasi-maximal slicing condition supplies an elliptic equation for a. Although elliptic equations can become 
computationally expensive, multigrid techniques, if applicable, can reduce significantly the cost. The same choice 
was made by Choptuik et al 5]. Like those authors we encountered several difficulties, the most severe being that 
in near-critical collapse situations the discretized version of the energy constraint equation for ip failed to remain 
diagonally dominant so that the multigrid algorithm failed to converge. Garfinkle et al ^} used a conjugate gradient 
method which gets round this particular problem. If instead we used an evolution equation for tp there appeared to 
be a serious violation of the energy constraint. 

We decided therefore to look at the other extreme, to minimize the number of elliptic equations solved. There are 
many, apparently different, ways of doing this, but we chose to implement the Z4 extension of Einstein's equation 
suggested by Bona et al @], Q, but adapted to the (2+l)+l scheme. The Z(2+l)+l scheme is outlined in section 
III1I There the energy and momentum constraints as well as the Geroch constraint (which arises in the (2+1) + 1 
decomposition) are converted into evolution equations for the new variables. 

Pursuing this philosophy we would like to impose dynamical gauge conditions, which are discussed in section llVl 
These are generalizations of the harmonic gauge condition. However cylindrical polar coordinates are not harmonic, 
and this means that we can only do this for the lapse function. For simplicity we assume that the shift vector vanishes. 

We now have a system of pure evolution equations and it is highly desirable to cast them into first order form. This 
is carried out in section W\ where we express them in conservation form with sources. We show in section IVT1 that this 
system is strongly hyperbolic and we write down the conversions between conserved and characteristic variables. The 
characteristic speeds (divided by the lapse function a) are 0, ±1 and iy/J where / is a constant parameter in the 
lapse evolution equation. (The choice / = 1 gives a harmonic time coordinate and a symmetric hyperbolic system.) 

Knowing the characteristic structure enables us to formulate constraint-preserving boundary conditions for the 
initial-boundary value problem in section IVIII 

In section lYlIII we return to our starting point, regularity on axis. By a judicious choice of new dependent variables 
(which does not affect the hyperbolicity) we can write our first order strongly hyperbolic system in a form where each 
and every term is manifestly regular on axis. 

Finally we discuss some of the implementation possibilities and applications of our scheme in section llXl 



II. THE (2+l)+l FORMALISM 

Spacetime is assumed to be a four-dimensional manifold {AA,g a p) with signature ( — h H — h) and a preferred polar 
coordinate chart (t, r, z, ip). Axisymmetry means that there is an everywhere spacelike Killing vector field £ = d/dip 
with closed orbits. The idea is to perform first a Kaluza-Klein-like reduction from AA to the Lorentzian three- 
dimensional manifold Af formed by the trajectories of the Killing vector. (This was carried out for vacuum spacetimes 
by Geroch . Nakamura et al [La. fig extended the reduction to include general matter sources and then considered 
the case of a perfect fluid. Choptuik et al [f| added a massless scalar field source, albeit without rotation. The current 
analysis includes arbitrary sources and rotation.) After this an ADM-like reduction is applied to Af. We shall follow 
the notation of Nakamura et al. [lH Il5l | , with some, clearly stated, changes. In the following, Greek indices range 
over t, r, z, if, lower-case Latin indices over t, r, z, and upper-case Latin indices over r, z. 

Tensor fields M""^.., in A4 which are both orthogonal to £ and are axisymmetric, i.e., 

M a -p„£0 = M- = . . . = 0, C<U"- = 

are said to be tensor fields in Af. Some basic tensor fields in Af are the norm of the Killing vector 

A 2 = 9a^ > , 

the (Lorentzian) metric in Af, 

h a /3 = gap — A~ 2 £ a £,9 , 
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the Levi-Civita tensor 



and the "twist vector" 



£a/37 — A 1 £ a (3'yS^ S j 



oJ a = e a ^V^ s , (2) 

where V is the covariant derivative of g a p. Note that neither — S a v nor = g Q¥ , is in A/". However if a solution 
of the equations in Af is given then the solution in M. can be reconstructed (Tj • We shall argue that a reconstruction 
is not necessary. All physically relevant quantities in M. have their counterparts in M . 
We define the projections of the energy-momentum tensor as follows, 

r a = \- 2 h a ^% s , (3) 
T a p — h a 1 hp S T 1 s. 

Note that near the axis A = 0(r), and the factors of A are included in © to ensure that the projections are either 
0(1) or 0(r) near the axis, so that regularity conditions (cf. section |VIII|) are easily imposed. 

Following the standard ADM procedure, the three-dimensional manifold N is foliated into two-dimensional spacelike 
hypersurfaces E(f) with unit normal n a — —ad a t and intrinsic metric 

H a b = hub + n a rib ■ 

The line element in J\f takes the form 

ds 2 = -a 2 dt 2 + H AB (dx A + r] A dt) (dx B + r] B dt) , 
where a is the lapse function and r/ A is the shift vector. The extrinsic curvature is given by 

Xab = —Ha C Hb d ri( c \d) = —jCnHab , 

where | denotes the covariant derivative compatible with h a b- The trace of the extrinsic curvature is abbreviated as 
X = XA A ■ Further variables defined in each S(t) are the alternating symbol 



CAB = Tl CcAB ■ 



V = -A"Va., 



the (/^-component of the extrinsic curvature, 



and the rotational degrees of freedom, 

E A = X - 3 e Ab u; b , 
= \- 3 n a uj a . 

Again, the last two definitions differ from those in by factors of A, and we write B v with an upstairs index. 

The various projections of the energy-momentum tensor are 

J* = -n a r a , 

S A = H a A r a , 

PH = n a n b T ab , (4) 

J A - -H A a n b T ab , 

Sab = H A a H B b T ab , (5) 

and, of course r defined in © • 
For convenience we also introduce 

Aa=oT X oc,a, L a = \~ x \ : a- (6) 
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We can now write down the projections of Einstein's equations 

which split into sets of elliptic constraint equations and hyperbolic evolution equations. 
The constraint equations are C = C A = C v = 0, where the energy constraint is 

C = i( X 2 - XABX AB + {2) R) - L A U - L A L A + - \X 2 {E A E A + B^ 2 ) - K p H , (7) 

with || denoting the covariant derivative of Hab, ^Rab the Ricci tensor of Hab and — ^Ra a the Ricci scalar. 
The linear momentum constraints are 

Ca = XA B \\ B - (x + K^) >A + L b xab - L A K V V - \X 2 B*e AB E B - n.J A , (8) 

and the angular momentum or Geroch constraint is 

C v = \E A \\ A + lL A E A - nJ^ . (9) 

We shall express all the evolution equations in terms of the Lie derivative along the normal direction, 

£„ = a^ 1 (d t - C n ). 

The time evolution of the metric variables is given by 

C n H AB = -2xab, (10) 
C n X = -XK^ . (11) 

The evolution equations for the extrinsic curvature variables are 

CnXAB = {2) Rab + (x + K ip v )xab - ^Xa° Xcb - A A \\ B ~ A A A B - L A \\ B - L A L B 
-\X 2 [e Ac e BD E c E D - H AB (E C E C - B* 2 )] 

-k [Sab + \H AB (p H - S c c - r)] , (12) 
£ n K v * = (X + K^) - L A ]]A - L A (L A + A A ) - \X 2 (E A E A - B* 2 ) 

-\k( P „-S a a +t) . (13) 
The rotation variables evolve according to 

C n E A = {x + ZK^)E A + e AB B^ tB +e AB (A B + 3L B )B^ -2kS a , (14) 
C n B* = X B v + e AB E A A B +e AB E Al]B . (15) 

(The principal parts of these equations resemble the axisymmetric Maxwell equations, which justifies the notation. 
The Geroch constraint @ is also a "Maxwell equation".) Energy-momentum conservation V a T a/3 = implies the 
following evolution equations for the matter variables: the energy conservation equation 

CnPH = -J A \\ A - J A {2A A + L A ) + X PH + XabS ab + K^r + X 2 E A S A , (16) 

the Euler equations 

C n J A = -S AB ^ B + J A (x + K^)~S AB {A B +L b )-A aPh + L a t 

+X 2 (E A J* + e AB S B B*) 1 (17) 
and angular momentum conservation 

C n J* = -S\ A - S A (A A + 3L A ) + J^( x + 3K V *) . (18) 

In many situations there may be a conserved particle number density N a satisfying C^N a — 0. If we introduce the 
projections 

v a = h a p N p , 
S A = H A b v\ 
a = -v a n a , 

then particle number conservation V a N a — implies the following evolution equation for a, 

C n a = - T, A (A A + L A ) + a( X + K v ") . (19) 
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III. THE Z4 EXTENSION 

Bona et al. suggested the following covariant extension of the Einstein field equations: 

R a p + 2V (a Z m = k (T a/} - \Tg aP ) , (20) 

which reduces to the Einstein equations if and only if Z a = 0. 

We shall assume that Z a shares the axisymmetry, i.e., C^Z a = 0. To apply the (2+l)+l formalism to equations 
(|2U[I . it is convenient to rewrite them as Einstein's equations G a p — nT a with a modified energy-momentum tensor 

~ 2 

T a = T a p (V(o,iTg) — \g a f3^ 'yZ 1 ) ■ (21) 

The Z covector is decomposed into parts parallel and orthogonal to the Killing vector £™, 

Z v = A 2 £ a Z a , Z a = h a a Z a , 
and further parallel and orthogonal to the timelike normal n a , 

9 = -n a Z a , Z A = H A a Z a . 
We now compute the modified matter variables corresponding to T a p, finding 

t = r + k- 1 [£ n 6 + Z A U + (A A - L A )Z A + [K v * - X )B] , 
Sa = Sa + k,- 1 [-Z* A + B*e AB Z B + E A 6] , 

J* = J V + K- 1 [CnZ* + E A Z A ] , 

Sab = S AB + ^ 1 [-2Z (Am + 2 X AB0 (22) 
+H AB {C n 6 + Z C \\ C + {A c + L C )Z C - (x + K V *)B} ] , 
J A = Ja + n- 1 [C n Z A - O.a + 2 XAB Z B + A A e] , 

" V 



Ph = PH + ^ 1 [C n 9-Z A u + (A A -L A )Z A + { x + K^ 



Inserting the above into the standard (2+l)+l equations, we obtain the Z(2+l)+l equations. Because of the Lie 
derivatives in (|22|) . the constraints (7H5|) become evolution equations for the Z covector, 

C n 9 = C + (L A - A A )Z A + Z A U - ( x + K v *)6 , (23) 
C n Z A = C A - 2 XAB Z B - A A 9 + 9,a , (24) 
C n Z v = C V -E A Z A . (25) 

The main evolution equations are modified as follows, 

£nXAB = ■ ■ ■ + 2Z ( A \\ B ) — 2xab0 , 
C n K v v = . . . + 2L A Z A - 2K,^9 , 
L n E A = ... + 2Z ip ' A -2E A 9~2B v e AB Z Bl 

where ... denote the right-hand-sides of iff^) . ifPS)) and (|T4"|) . respectively. The remaining evolution equations are 
unchanged. 
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IV. DYNAMICAL GAUGE CONDITIONS 

To complete our evolution formalism, we need to prescribe the gauge variables a and r\ A . It is well-known 0,0 
that the Einstein equations can be written in symmetric hyperbolic form 

5 75 5a/3, 7 ,5 ^ 

(~ denoting equality to principal parts) by adopting harmonic gauge, which for the Z4 extension of the field equations 
becomes 

g-i 5 x a ]l5 = 2Z a , 

where we are treating x a as a scalar field. It would therefore be desirable to choose this or a related gauge for 
the system at hand. However, these gauge conditions are incompatible with cylindrical polar coordinates; the wave 
operator contains a term r~ 1 d r , and the right-hand-side of the evolution equation for the shift vector becomes singular 
on the axis. However, we can keep the time component of the harmonic gauge equation, which in (2+l)+l language 
reads 

d t a = -a 2 { X + K^ -29). 

This condition has been generalized to Q 

d t a = -fa 2 { X + K^-mO), (26) 

where / and m are free constant parameters. 

We choose the shift vector to vanish: T) — 0. Apart from simplifying the evolution equations, this choice also 
makes it much easier to set up outer boundary conditions because a variable shift would mean that the characteristic 
speeds could change sign at the boundary (see section IvTjl . Our choice for the shift vector means that henceforth 

C n = a^dt ■ 



V. FIRST-ORDER REDUCTION 



The Z(2+l) + l equations are first-order in time but second-order in space. In order to analyze their hyperbolicity 
and to be able to treat them with standard numerical methods for hyperbolic conservation laws, we perform a 
first-order reduction. We focus on the geometry evolution equations here; for given geometry, the matter evolution 
equations form a separate strongly hyperbolic system (see appendix lD|) . 

To eliminate the second spatial derivatives, we introduce new variables for the first spatial derivatives of the two- 
metric, 

Dabc = \QaHbc i 

and we regard A a and La defined by JSJ as independent variables. Evolution equations for these new first-order 
variables can be obtained from 11U|) . and (|26|l by commuting space and time derivatives: 

£ n DABC — ~Xbc,a , (27) 
C n L A = -i<W (28) 
C n A A = -f{ X ,A + K^ A -me }A ). (29) 

The two independent traces of Dabc are denoted by 

D 1 a = Dab B , D 11 a = D B ba , 

where indices are (formally, for D is not a tensor) raised and lowered with the 2-metric Hab ■ The De Donder-Fock 
decomposition 0,13 is used for the Ricci tensor, 

^Rab = —D c ab.c + 2D 11 \a,b) — D 1 (A,B) — ^D C abD iic 

-T cab (2D iic -D ic )+ AD C daD CD b - T A cdT b CD , 



where of course the Christoffel symbols are given by 

Fabc = D C ab + D B c a - Dabc- 
It is then straightforward to write the Z(2+l)+l equations in conservation form 

8 t u+ [aT D {u)] D =aS(u). 

Here, the set of conserved variables is 

u = (Hab, A, a, Dabc, L a , A a , xab, K^,E a ,B*, 9, Z A , Z*) T , 
and the components of the flux vector are given by 



J~ B AB 




o, 




T D 

■r a 




o, 




J~ a 




o, 




*~ Dabc 




SaXbc , 




tD 
J- L A 




5aK^ , 




J~ A A 




Sa fix + KJ> - wff) , 




xab 




D D ab - Sf A (2D H B) 


+ 2Z B) 






L D , 




J- E A 




-2H AD Z V - e AD B v . 




J- Bv 




-e AD E A , 




J- 




D ID _ D IID +L D _ 


z D , 


J- Z A 




-XA D + S°( X + K^ 


-o), 


tD 

J~ Z'f 









D 1 B) - L B ) - A-i 



5(u) is a source term containing no derivatives, 
Sbab = -%Xab , 

S a = -fa( X + K v *-m9), 

Sdabc = °) 

Sl a = 0, 
Sa a = 0, 

S XAB = AcD c ab+A (a (-2D ii B) +D i b) +L b) +A b) -2Z b) ) 

-L A L B - A A A B + T C ab(L c + A c ) + xab(x + K^>) - 2 X a C Xcb 
-\\ 2 [e A ce B DE c E D - H AB {E C E C - B^)} 
-k [S AB + \H ab (ph - S c c - r)] - Y CAB (2Z C + 2D IIC - D IC ) 
—29xab — 2D C abD h ° + 4D CDA D CD B — ^ acd^ b CD , 
S Kv * = K^( X + K^)+La(2Z a -L a -D ia )-\\ 2 {E a E a + B^ 2 )-2K^6 
-\k(ph - S c c + t) , 
S e a = (x + 3K v v -26)E A + e AB B v {iL B -2Z B + D 1 B ) - 2kS a + (4D iia - 2A J 

S b , = X B V + e AB E A D I B , 
S e = A a {D ia -D iia + L a -Z a ) + DabcD abc -\T A bcT ABC -\D J A D IA 

-L A (L A + D ia ) + i( X 2 - xabX AB ) + XK^ - \\\E A E A + B^) - np H 

+ (L A -A a + D i a )Z a - (x + K v *)0 , 

Sz A = -A bX a B +A A (x + K^ -6)+xab(D ib +L B -2Z b )-T c abXc B 

-L A K v f - \\ 2 B^ AB E B - kJ a - A A 9 , 
S z * = ±E A (D I A + 3L A -2Z A - A A )- kJ* . 
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Note that Hab, A and a have vanishing fluxes and thus trivially propagate along the normal direction. The rotation 
variables (E A ,B (p ,Z (p ) form a decoupled subsystem analogous to Maxwell's equations. 

VI. HYPERBOLICITY 

To investigate the hyperbolicity of the Z(2+l)+l system, we pick a unit covector [ia and define an orthogonal 
covector 

tta = £abH B , (32) 

so that 

pL A H A = TTATr A = 1 , \IaK A = . 
Thus (h a ,tt a ) form an orthonormal basis for T(S). Projection along fi and tt is denoted as 

y 11 ee V a ha , v 1 - = V a -k a ■ (33) 
It can be verified that the Jacobian matrix 



du 



corresponding to the flux in the /i-direction is real, diagonalizable and has complete sets of left and right eigenvectors 
for arbitrary ha, i-e., the system is strongly hyperbolic. The characteristic speeds are 3 

A = 0, 

Xf = ±1 (speed of light) , 
A/ = ±V7 (gauge speed) . 

For / ^ 1, all characteristic speeds are causal, for / = 1 (harmonic slicing) they are all physical. 
The characteristic variables (left eigenvectors dotted into the conserved variables) are given by 

• Normal modes (eigenvalue 0) 

I 
I 
I 



'o.i 


= fm(D llxx 






~ z \\) 


~f(D - 




+ L|,) + A||, 


'0,2 


= fm(D xu 


+ L ± 


- D \\\\x 


-Zx) 


-f( D M\\\ 


+ D XX , 


i + L X ) + A x 


'o.3 


= D M\h 














'o.4 


= D XX|| , 














■0,5 


= D ±±± , 














'n.6 


= L x , 














■0,7 


= A ± , 















along with the zeroth-order variables Hab , A and a 
• Light cone modes (eigenvalue ±1) 



if tl = e±(D ll±x + L ll -D xxll -z ll ), 

lf a = X\\±±\{Ax + D m -D xxx + L x -2Z x ), 

h,S = X-L-L ± -D||_LX > 



n,4 — Xl *> 



3 Note that a factor of a has been taken out of the fluxes 1311 . Note also that for a nonzero shift vector r/ A , a term —a 1 r/H would have 
to be added to the above characteristic speeds. 
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Lapse cone modes (eigenvalue ±\/7) 



If = A \\ ~ 1 ^H-L-L + L W - D ±M\ - Z W ' 

7("»-2) 



±V7 



Xiiu + Xx± + K v v -[ J J_ -' + 2)9 



(34) 



The inverse transformation from characteristic to conserved variables is given in Appendix iBl 

In the special case / = 1, we must set m = 2 to enforce strong hyperbolicity, and the singular term (m — 2)/(/ — 1) 
in (|34H is to be replaced with an arbitrary fixed constant (e.g. for simplicity). In this case, the system is even 
symmetric hyperbolic, i.e. 2?" is symmetrizable with a symmetrizer that is independent of the direction fi A . This is 
not surprising because the choice of parameters / = 1, m = 2 in (|26[) corresponds to harmonic slicing. An energy for 
the system is given by 

£ = xabX AB + ^abc\ ABC + (K^+x-20) 2 +A a A a 
+V A V A + K^ 2 + L A L A + E A E A + B v2 + AZ^ 2 , 

where 

V A = A A + D* A + L A - 2D 11 A - 2Z A , 
A A bc = D A bc + 5 a b Vc)- 

£ is positive definite and its time derivative is (considering principal parts only) a total divergence. 

VII. CONSTRAINT-PRESERVING BOUNDARY CONDITIONS 

In the Z(2+l)+l formalism, the standard (2+l)+l constraints (00 are replaced with the algebraic constraints 

6 = Z A = Z v = . (35) 

If those hold at all times, the (2+l)+l constraints are automatically satisfied by virtue of the evolution equations for 
the Z vector Q23H25fl . On the initial time level, one imposes l|35|l and solves the (2+1) + 1 constraints 00 so that 
both the Z vector and its time derivative vanish. In addition, the outer boundary conditions have to be chosen so 
that no constraint-violating modes propagate into the computational domain at any time. 

A propagation system for l|35|) follows from the contracted Bianchi identities, which imply V a T al3 — 0. If we 
impose the standard matter evolution equations V ' a T a @ = in addition, l|21|) yields the following homogeneous wave 
equation for the Z covector: 

VpV p Z a + R a0 Z p = . (36) 

To obtain the (2+1)+ 1 reduction of (|36|l . we take an additional time derivative of the evolution equations for the 
Z vector (|23II25|I . To principal parts, we have 

C?J-d B d B e ~ 0, 
C 2 n Z A -d B d B Z A ~ \d B [8 A (A B +D IB + L B )-d B (A A + D IA + L A ) 

-2(d c D BC A - d B D n A )] , (37) 

£ 2 n Zf-d B d B Z v ~ 0. 

Note that the right-hand-sides are zero if the ordering constraints 

9dD A bc = d A D DB c , 

d B L A = d A L B , (38) 
d B A A = d A A B 

are satisfied, which corresponds to commuting of partial derivatives in the definition of the first-order variables. 
Although this holds analytically, it may not obtain numerically, and so the terms have to be included in the subsidiary 
system |[37J. 
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We choose an orthonormal basis fi A , ir A such that [i A is normal to the outer boundary under consideration. Indices 
contracted with \i A (n A ) are denoted by |j (_L). Equations 1)3 7|) decompose into 

c 2 j - afe - d 2 x e ~ o, 

C 2 n Z l{ - 3\Z\\ - d\Z\\ ~ d\D^ 

+±8 {l d ± (A ± - D ±n + D ±±± - 2D n± + L±) 

+ld 2 ± (-A ]l +D - - , 

C 2 n Z ± - 8\ Z ± - d 2 ± Z ± ~ ±d*(-A ± -D m +D ±±± -L ± ) 

+±d {l d ± (A {l + D - D ]l±± - 2£»xx|| + L {{ ) 

,2 



+d±D U ±, 
C 2 n Z v - d^Z^ - dlZ v ~ 0. 

This second-order system can be written in first-order form by introducing new variables for the temporal and spatial 
derivatives of the fundamental variables. Absorbing boundary conditions can then be imposed in the || direction, 
which read (after replacing the new variables with their definitions as first derivatives of the fundamental variables) 

C n 9 = -d l{ 6, (39) 

C n Z {l = -d ll (Z ll +D ±±l[ )-±dx(A ± -D xn +D ±±± -2D u± +L ± ), (40) 
C n Z ± = -d\\ [Z±-^A x + D ±m -D ±X j_ + L ± )\ 

+D - D [l±x - 2D ±±ll + L||) , (41) 

£ n Z* = -d {{ ZV, (42) 

where = denotes equality at the boundary. 

The straightforward way of setting up stable outer boundary conditions for the main evolution system is to set all 
ingoing modes to zero at the outer boundaries of the computational domain (so-called absorbing boundary conditions). 
For a symmetric hyperbolic system, this leads to a well-posed initial boundary value problem (IBVP) [rH lisj . Our 
goal is to investigate under which circumstances such boundary conditions imply the constraint-preserving first-order 
conditions I|39II42|I . Suppose l~ is an incoming mode with characteristic speed —A (A > 0). Setting l~ to zero at the 
entire boundary at all times implies that 

o = d x r = Li(r) (43) 

and 

o = £ n r ~ ASyr-ax^V =L 2 (r), (44) 

where = denotes equality at the boundary and we have used the general general evolution equation for l~ . The 
relations (|43|l and (|44|l will now be used to manipulate the general evolution equations for the Z vector at the 
boundary. As mentioned in section IVTl the evolution system decouples into a non-rotational and a rotational part at 
the level of principal parts, and we consider the two subsystems separately. 

Let us start with the general evolution equations for 9 and Za, decomposed with respect to our boundary- adapted 
basis: 

C n 9 ~ -d||CD||xx-£>xx|| +L\\ -Z||)-dxCDx|||| -D u± + L±-Z ± ), 



C n Z ± ~ d\\X\\± - d±(x\\\\ + K< 



If we add to these suitable combinations of equations (|43J) and 144(1 , 

c n e += LaHlx), 

£-n,Z\\ += L2{ — l ll +l lz + l lA ) + Li( — l l ^), 

C n Z±_ += L 2 {-1 12 ) + Li(li tl - li 3 - IJ ) , 

we obtain precisely the constraint-preserving boundary conditions (|39H41|I . provided that we choose harmonic gauge 
f = 1, m = 2, — 0. (For different values of / and m, the derivatives tangential to the boundary in the evolution 
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equation for Z±_ differ from the constraint-preserving boundary condition Q41[l.l Hence absorbing boundary conditions 
for the non-rotational part of the evolution system preserve the constraints. 

However, this is not true for the rotational subsystem. A little calculation shows that the only dissipative boundary 
condition consistent with (|42[l is 

Fortunately, there is a much simpler solution to this problem: we delete Z v entirely from the evolution system. Note 
that setting Z v = does not break general covariance because our spacetime has a Killing vector and we simply 
choose Z a t; a = 0. The remaining rotational subsystem is still symmetric hyperbolic, and the constraint 

Lstp — 2 rj ,A 

has zero speed with respect to the boundary: 
The absorbing boundary condition 

= E 1 - + B v = (45) 

leads to a well-posed IBVP for this system. 

In addition to the Z constraints, there are differential constraints related to the definition of the first-order variables, 

Dabc — \9aHbc i 
L A = A-^aA, 
Aa = a~ x dAa, 

and the ordering constraints (|38|l . By virtue of the evolution equations (|10l 1111 l2l))l and (|27l 1281 l2T)|) . it is easy to see 
that the differential and ordering constraints have zero speed with respect to the boundary. 

In summary, we have shown that for harmonic gauge, absorbing boundary conditions for the non-rotational part 
of the main evolution system are automatically constraint-preserving, and that there are two ways of setting up 
constraint-preserving stable outer boundary conditions for the rotational subsystem. We thus have a well-posed 
IBVP with constraint-preserving boundary conditions. The result concerning the non-rotational part carries over 
immediately to the general Z4 system without spacetime symmetries. Constraint-preserving boundary conditions for 
the Z4 system have recently also been analysed by Bona et al. Q: in essence, those authors propose to implement 
the constraint-preserving boundary conditions (|39H41|) directly. While it is not clear analytically that such boundary 
conditions lead to a well-posed IBVP, Bona et al. are able to prove a necessary condition and to demonstrate numerical 
stability. 



VIII. REGULARITY ON AXIS 



Extra care is needed in axisymmetric situations because the coordinate system that is adapted to the symmetry, 
cylindrical polar coordinates, is singular on the axis. This leads to rather strong regularity conditions on tensor fields, 
as explained in appendix ^ In addition to the axisymmetry, we would like to impose reflection symmetry about the 
z-axis so that we only need to evolve one half of the (r, z) plane. This implies that tensor components are either odd 
or even functions of z. 

Let us first deal with one of the regularity conditions for 2-tensors M a p, which follows from 

4fe^ = l + 0(r 2 ) 



near the r-axis. For the metric g a [3 this implies 

A 2 g w 



1 + 0(0 



We enforce this condition by replacing A with a new variable s defined by 
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To satisfy the corresponding condition for the extrinsic curvature, we introduce Y via 



H„ 



r 2 Y 



(47) 



(note that K, 



Vl p — \ 2 K ip ip ). Similary for the energy-momentum tensor, we set 

t = h r t . 



We remark that the definitions of the variables s l|4f)|l and Y 147(1 can be viewed as a generalization of those in Q and 

The remaining dependent variables are redefined by taking out systematically the leading order of r and z: 

H rr H rr , H rz — H rz j (rz) , H zz H zz , 

D rrr — 2 rr,r j ** , D rrz — 2 rz ' r /^*' J~) rzz — 2 zz - r ' 

f) — -H Iz D —-H Iz D =-H Iz 

zrr — 2 rr,z / ~ 3 zrz 2 rz,z / ^ , ^ zzz 2 zz -. z ' 1 

S r s r /r , S z S _ z I Z , 

a = a , A r — a^ 1 a^ r /r, A z —a^ 1 a, z /z, 

Xrr = Xrr, Xrz = Xrz I \rz) , \zz = X zz , 

E r = E r /r, E z =E z /z, B v =B v /(rz), 
6 = 6, Z r = Z r /r, Z z = Z z /z, Z* = Z* , 
a = a , p K = p H - a , J r = J r /r , J z = J z /z , J v = J v , 
E r = E r /r , Yi z = S z /z , S r = S r /r , S z — S z / z , 
S rr S rr , S rz — S rz f (tz*) , S zz — S zz . 

It can now be verified 4 that the Z(2+l)+l equations can be written in terms of the new variables u in the form 



dtU 



aT {r } (u) 



aS(u) , 



(48) 



where the fluxes T and the source S are even in r and z and manifestly regular on the axes (cf. appendix IC| . 

The above transformation of variables does not affect the hyperbolicity of the system. To see this one should note 
that 



D r 
D, 



- D rrr f t , D rrz 

D zrr I % •> D zrz 

D 



D rrz /{r 2 z) , D 

rzz 

D zrz /(rz 2 ) , D zzz 
D zr 



Xrz 



A r ~ A r I r , A. 
Xrz/ (rz) , x*z - Xzz, Y ~ [K^ - 
E r 



Drzz/r, 

D zzz /z, 
' A z /z, 

~r \ I 2 

■)/r , 



H, 



E r /r, E z ~E z /z, B v ^B v /(rz), 
8-8, Z r ~Z r /r, Z z ~Z z /z, Z v ~ Z v , 

Jr — Jv/t , Jz — Jz/z, 

S r ~ S r /r , S z ~ H z j z , S r — S r / r , S z — S z j ' z , 

S rr — S rr , S rz — Srz/ (yz) , S zz — S zz , 

where ~ denotes equality to leading derivative order (the lower-order terms are absorbed in the source terms) 5 . Hence 
the variables occuring in the fluxes undergo a linear transformation, which does not affect the hyperbolicity and leaves 
the characteristic variables unchanged. 



4 The calculations are rather lengthy. We used the computer algebra system REDUCE. 

5 Note that H rr has zero flux and is thus to be treated as a background quantity. 
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By writing the equations in terms of derivatives with respect to r 2 and z 2 instead of r and z, however, we have 
ultimately to multiply the fluxes by 2r and 2z, respectively, because 

d r 2 — — d r 
2r 

and similarly for z. This means that in (r 2 ,z 2 ) coordinates, the characteristic speeds behave like ~ r and ~ z, 
respectively. For computational purposes it is highly desirable to avoid the non-uniform characteristic speeds which 
would occur on a (r 2 ,z 2 ) grid. We therefore recommend working in the original (r, z) grid and discretizing the 
derivatives in (|48|l with respect to r 2 and z 2 . This has another, perhaps more fundamental, advantage: in (r, z) 
coordinates, we can enforce Neumann boundary conditions for all our variables on the axes (since they are even 
functions of r and z), whereas in (r 2 ,z 2 ) coordinates the boundary conditions on the axes are not known. The 
authors of faced a similar problem but decided to work on an (r 2 , z 2 ) grid, using extrapolation on the axis, which 
led to numerical instabilities. 

Choptuik et al p| and Garfinkle et al take a slightly different approach. They include terms which are formally 
irregular on axis, but which are regular once appropriate boundary conditions are imposed. Following an idea of 
Evans Q these terms are differenced with respect to r 2 or r 3 . 

The transformation from conserved to primitive fluid variables and vice versa can be made manifestly regular on 
the axes by redefining the velocities 

v r = v r /r, v z = v z /z 

and leaving the remaining primitive variables unchanged. Note in particular that the modified general matter variable 



f = r- 2 ( t -^l)= P hW 2 (e 2 ^H rr v* 2 - (49) 

\ tl rT J \ £l rr J 

is then automatically regular. 

To compute the characteristic variables (geometry and matter), one can start from the regularized conserved 
variables u, compute the original conserved variables u and evaluate the characteristic variables (section IVII and 
appendix ID 3(1 . While this transformation is perfectly regular on the axes, the inverse transformation involves negative 
powers of r and z. In order to implement outer boundary conditions at r = r rnax (z — z max ), one has to ensure 
that the inverse transformation in the r-direction (z-direction) is regular at z = (r = 0). This can be achieved by 
redefining the characteristic variables in a similar fashion as was done for the conserved variables above. We choose 
to make the following replacements 

• for the characteristic variables in the r-direction: 

^0,4 ~ > ^0,4 = Z 2 ^0,4 + ^J~fl ^°' 5 ^ ' 

and the leading order of z is taken out of the remaining characteristic variables 

• for the characteristic variables in the z-direction: 

; T -2 ( 1 , r ^rz , \ 

'0,4 — * <0,4 — r I '0,4 + rj^ to, 5 I i 



^0,6 - * ^0,6 — r 3 (^0,6 — '0,5) , 
'1,4 1,4 — ' V-1,4 b l,3) ' 

and the leading order of r is taken out of the remaining characteristic variables 

We have checked with REDUCE that the transformation from characteristic variables to conserved variables and vice 
versa is then manifestly regular at the entire outer boundary. 

However, the transformation from characteristic variables in the r-direction (z-direction) to conserved variables is 
still formally singular at r = [z = 0), and no regularization procedure can cure this problem. To understand this, 
one can note that the characteristic variables in the r-direction (z-direction) do not have a definite parity in r (z). For 
this reason, numerical methods operating in the space of characteristic variables (typically ones based on the solution 
of the Riemann problem) appear to be unusable near the axes. 
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IX. CONCLUSIONS 

We started out by trying to understand why so many previous attempts to evolve axisymmetric spacetimes failed 
because of on-axis instabilities. This led to a detailed study of the behaviour of the components of axisymmetric 
tensors with respect to cylindrical polar coordinates, given that the components with respect to cartesian coordinates 
were regular in the neighbourhood of the axis. This s ugg ested, most strongly, that we should make a Kaluza-Klein-likc 
reduction to a 3-dimensional Lorentzian spacetime |llj . a step taken earlier by Maeda et al 0], [l^. (A different 
version of this idea, without rotation, has been pursued by 0|.) 

Choptuik et al's evolution [j| was fully constrained — elliptic equations were solved for the metric components, 
making use of multigrid for computational efficiency. We made similar choices and obtained similar results. However 
in strong field situations, e.g., near-critical collapse of Brill waves, the diagonal dominance of the matrix describing 
the discretized elliptic operators failed, and so the multigrid iterations failed to converge. 

We have not found a satisfactory resolution of this problem, and so we went to the other extreme — our evolution 
algorithm involves no elliptic operators. To achieve this we have used the Z4 formalism 0, adapted to our axisym- 
metric reduction, an evolution equation for the lapse function and a zero shift gauge. Our first order hyperbolic system 
of conservation laws (with sources) is strongly hyperbolic and, for one choice of parameters, symmetric hyperbolic. 
In addition the dependent variables have been carefully chosen so that each and every term in the system is regular 
on axis. We have allowed for arbitrary matter sources and have presented explicitly the details for a perfect fluid. 

Initial conditions have to be imposed via the constraint conditions, which means solving elliptic equations. By 
means of conformal rescalings on the initial hypersurface one can avoid the diagonal dominance problem on that 
hypersurface and so multigrid techniques are applicable. 

One must choose a high resolution shock capturing method for the evolution of the hyperbolic system. Our system 
contains one subtle disadvantage if we choose to work with numerical algorithms that require a transformation between 
conserved and characteristic variables. The conversion from characteristic variables to our "regularized" conserved 
variables, see section lYUH is formally singular on the axis. However, we have shown that the transformation can be 
used to impose outer boundary conditions in a regular way. 

Because the evolution does not involve elliptic operators it is straightforward to introduce adaptive mesh refinement, 
where appropriate, to enhance the resolution. 

We shall be using our formalism to evolve a number of interesting axisymmetric spacetimes, including the effects 
of rotation and perfect fluid sources. 
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APPENDIX A: CONSEQUENCES OF AXISYMMETRY 

We want to use a (t,z,r,ip) chart adapted to the Killing vector £ = d/dip. Unfortunately this chart is singular at 
the axis r = 0. We shall assume elementary flatness: in a neighbourhood of the axis we can introduce local cartesian 
coordinates x A = (x,y) such that 



x 



x = r cos tp, y = r sin tp r — \J x 2 + y 2 , tp = arctan ■ 

With respect to cartesian coordinates the Killing vector is 

£ = —yd/dx + xd/dy. 

Notice that this representation is valid everywhere, while £ = d/dip is valid only for r > 0. 

We say that a scalar function fix ) is regular on axis if / has a Taylor expansion with respect to x and y about 
x A = convergent in some neighbourhood of r = 0. (Throughout this section we are ignoring t, z dependencies which 
are implicit in all calculations.) If / is axisymmetric 

C 6 f = = df/dtp f = f(r), r>0. (Al) 

If / had a Taylor expansion in r and was regular on axis then the expansion could contain only even powers of r 
since r — x 2 + y 2 , and r has no such expansion. While the conclusion is correct, the argument is fallacious because 
equation ilAlfl is invalid at r — 0. Instead we start from 

K = -yf. x + xf. v = 0, (A2) 



which is valid everywhere. In particular we may differentiate <|A2|1 an arbitrary number of times with respect to x 
and also with respect to y. Then setting x A = we obtain linear recurrence relations between the Taylor coefficients 
on axis. These can be solved to show / = J2 n fn{x 2 + y 2 ) n . 

Next consider a vector field u a . For a = (t, z), C^u a = implies du a /dip = 0. This reduces to the scalar field case 
and we may deduce u a = u a {r 2 ). For u x and u v we have 

du x /dp> + u y = 0, du y /dip - u x = 0. 

The general solution for r > is 

u x = a{r) cos ip — b{r) sin tp, u v = a(r) sin tp + b(r) cos tp. 
However in the cartesian chart we have 

~yu x , x + xu x /y + u v = 0, ~yu v iX + xu v iV - u x = 0. (A3) 
Clearly u A = on axis. We write a — ra etc so that 

u x — xa — yb 1 u v = ya + xb. (A4) 



We now regard a and b as unknown functions of x and y to be determined by substituting (|A4|) into (|A3|) , differentiating 
the latter an arbitrary number of times, and then solving the recurrence relations for the Taylor coefficients of a and 
b. Again we find that a and b are even functions of r. Thus in the (t, z, r, tp) chart an axisymmetric vector field which 
is regular on axis must take the form 

u a = (A,B,rC,D), (A5) 

where A, B, C and D are functions of t, z and r 2 . 

Next consider an axisymmetric covector field Lu a which is regular on axis. For a — (t,z), C^Lu a — implies 
duJa/dtp — 0. This reduces to the scalar field case and we may deduce ui a = oj a (r 2 ). For the other indices we find 

~yUx,x + %Ux,y +UJy = 0, -yCJy,x + XLUy^y -lj x = 0, (A6) 

which is equivalent to l|A3() . interchanging u A and to a- We therefore deduce the analogue of l|A4(l and hence, in polar 
coordinates, 

uj a = {E,F,rG,r 2 H), (A7) 



16 



where E, F, G and H are functions of t, z and r 2 . 

Finally we consider a symmetric valence 2 axisymmetric tensor field M a p which is both axisymmetric and regular 
on axis. For (a, b) = (t, z) we have C^M a b = and so M ao = M ao (r 2 ). Letting a — (t, z) we have 



-yM aX}X + XM aX: y + M ay = 0, -llM a y iX + X M ' ay >y ~ M aX = 0. 



(A8) 



This is essentially the same as <|A6ll and we may deduce M ar = rA a (r 2 ) and M av = r 2 B a (r 2 ). The remaining Killing 
equations are 

-yM xx%x + xM XX}V + 2M xy = 0, -yM yy , x + xM yy ^ - 2M xy = 0, -yM xy ^ x + xM xy ^ y + M yy - M xx = 0. 
If we introduce new variables u = \(M XX + M yy ), v = \{M XX — M yy ) and w = M xy then 

~V u ,x + xu -y — 
which implies u — u(r 2 ). The remaining equations are 

—yv iX + xv_ y + 2w — 0, —yw -x + xw -y — 2v = 0. (A9) 

For r > these can be written as 

v iV + 2w = 0, w ;tp — 2w = 0, 

so that 

v = a(r) cos 2tp — b(r) sin 2ip 1 w = a(r) sin 2ip + b{r) cos 2tp, 

where a and b are arbitrary functions of r. If we set a — a/r 2 and b — b/r 2 then 

v = (x 2 - y 2 )a ~~ 2xyb, w = 2xya + (x 2 - y 2 )b. (A10) 

(Note that (|A9(I and its first derivatives imply that v, w and their first derivatives vanish on axis, consistent with 
(|AT0|) .) Substituting (jA"T0|) into $9% gives 



x 3 a, y - x 2 y(a^ x + 2b, y ) - xy 2 (a, y - 2& x ) + y 3 a tX = 0, 
x 3 b iV + x 2 y(-b tX + 2a y ) + xy 2 (-2a tX - b >y ) + y 3 b tX = 0. 



(All) 
(A12) 



Differentiating these many times and, proceeding as in the scalar and vector cases we conclude that a and b are 
functions of r 2 . Thus 



Air 



u + (x 2 — y 2 )a — 2xyb, M yy = u — (x 2 — y 2 )a + 2xyb, M xy = 2xya + (x 2 — y 2 )b. 



Re-expressing these as polar components we obtain 



M„ = u + r 2 a, M rip = r 3 b, M w = r 2 (u - r 2 a). 



Finally combining all of the results we have 



M, 



'la/3 



I A B rD 

B C rE r 2 G 
rD rE H + r 2 J r 3 K 

\ r 2 F r 2 G r 3 K r 2 (H - r 2 j) ) 



(A13) 



(A14) 



where A,B,...,K are functions of t, z and r 2 
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APPENDIX B: TRANSFORMATION FROM CHARACTERISTIC TO CONSERVED VARIABLES 

The conserved variables are given in terms of the characteristic variables by 



D 


= -7/0,1 + ( 

+ Tf 0/ + 
ZJ \ 3 




= —-j, — lo,2 + 

jm 


D [l±x 


— 2V'l,3 '1,3 


Dx U 


= ^0,3 ; 


D ±M\ 


= ^0,4 i 


D ±±± 


= ^0,5 ; 


L \\ 


= a(^4 — 'l,4 




= fo,6 j 



( m ~ 2 ^/ 4-/ -4-7 ^ (/ m ~ 2 ) / ,1(1+ 1- \ 
(.'0,3 + '0,5 + k)fi) TT? '0,7 + 2\ l l 2 — 'l 2) ' 



2m 



2/m 



A {{ = 
A± = 
X\\\\ = 



f(m-2) n+ 
2(/-l) ( ^ 

fo,7 , 

/(m-2) 



2(/-l) 



+ 1 ) (£l + - 2^3 + ^3 + ^4 + ^4) 





+ 2V7 




X||_l = 




^,2) ' 


x±± = 


Wh + 


^1,3) ' 


= 




'1,4) > 


eH = 




'i; 6 ). 


£ x = 


1(^6 + 


^l,6J ' 




2 ^'1,6 


" llfi) . 


e = 


3(*£i + 


^l) , 


z \\ = 


- '0,4 + 


2 ( — ^1,1 + 'l,l + '£3 — 'l,3 + ^4 — 


z± = 


2"( ; 0,3 - 


Iq,5 + fo,6 + ^0,7) — 2(^1,2 — ^1,2) ! 


z v = 







Tensor components can then be computed as, for instance, 

Xab = Ha^bXWW + 2/u.( j4 7r_ B )X||-L + ^a^bX±-L 



APPENDIX C: REGULAR CONSERVATION FORM 



In this appendix we provide some details of the regular form of the Z(2+l)+l equations outlined in section IVIIII 
The equations are written in conservation form 1481) 



a^ 2) (u)l =aS(u), 



where the modified variables u are defined in section IVIIII 

The fluxes f {r \ ' and sources S have been obtained with the computer algebra system REDUCE and have 
been typeset automatically using the TeX REDUCE Interface (TRI). From the same source we generated C++ code 
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using the REDUCE source code optimization packages SCOPE and GENTRAN. The resulting expressions are rather 
lengthy. To get some idea of what is typically involved we state here the fluxes and source for the variable Y. 
Let us introduce the shorthand 

H = H rr H zz - r 2 z 2 H 2 rz 
for the determinant of the 2-metric. The fluxes are given by 

f ( p = 2H- 1 D rrr r 2 z 2 H- 2 H 2 z - 2H~ x b rzz - AH' 1 b zrr z 2 H^ H rz + AH^D^z 2 

+2H- 1 (-A r + 2Z r ) + 2H- 1 r 2 z 2 H rz (r 2 H- 1 H rz S r - s z ) , 
T ( p = 2H- 1 z 2 (-r 2 H rz S r + H rr S z ), 

and the source term is 

Sy = -kt + 2H- 1 x 2 rr z 2 H- 2 H 2 z + 2x rr H- r 1 Y + H^Xrv^H-J-H^H^Y - A Xr z) 
+H- 1 XzzH rr Y - AH^D^yz 2 !!- 3 !!^ 
-H- 2 D 2 rrr r 4 z 4 H- 3 H? z + ^R- x T) rrr t) rrz r 2 z 2 H~ 2 H rz 
+2H 2 b rrr b rrz r 4 z 4 H r 2 H 3 Z + AH l D rrr D zrr z 2 H r 2 H rz 
~\~2H D rrr D zrr v z H rr H rz 2H D rrT D zzz z H rz 

+H- 1 D rrr H 7 : 2 H rz (-5r 4 z 2 H rz § r - r 2 z 2 A r H rz + Ar 2 z 2 H rr s z - 2r 2 z 2 H rz s 

+Ar 2 z 2 H rz Z r + 2z 2 A z H rr - Az 2 H rr Z z + 2H rr ) 
+H- 2 D rrr r 2 z 2 H- 2 H 3 z (-r 4 z 2 H rz S r + r 2 z 2 H rr S z + 2r 2 z 2 H rz s + 2z 2 H rz + 2H, 



rr } 



-2H- 2 b rrz b rzz r 2 z 2 H rz - AH- x b vrz b zrr z 2 H^ - m- 2 b rrz b zrv r 2 z 4 k-}H 2 rz 

-\~Aiii D rrz D zrz v z H rz ~\~ D rrz D zzz z H rr 
+2H- 1 D rrz (2r i z 2 H- r 1 H rz S r + r 2 z 2 A r H^H rz - 2r 2 z 2 H~ r 1 H rz Z r 

-r 2 z 2 S z - z 2 A z + 2z 2 Z z - 1) 
+2H- 2 b rrz r 2 z 2 H 2 z (r' i z 2 H- 1 H rz S r - 2r 2 z 2 H^ 1 H rz ~s - r 2 z 2 s z - z 2 H^H rz - 1) 
+H 2 b 2 zz H rr + AH 2 b rzz b zrr z 2 H rz — AH 2 b rzz b zrz z 2 H rr 
+H- 1 b rzz (-r 2 S r -A r - 2.?) + H- 2 b rzz r 2 z 2 H rz {-r 2 H rz s r + H rr ~s z + 2H rz ~s) 
+2H- 1 b zrr z 2 (r 2 H- r 1 H rz S r - A r H~ x H rz - ~s z ) 
+H- 2 b zrr z 4 H 2 z (r 4 H- 1 H rz S r - 2r 2 H^H rz ~s - r 2 ~s z - AH„ l H rz ) 
+2H~ 1 b zrz z 2 (A r + 2s) + 2H~ 2 b zrz z 4 H rz (-r 4 H rz s r + r 2 H rr s z + 2r 2 H rz s + 2H rz ) 

+H- 2 b zzz z 2 H rr {r 2 H rz ~s r - H rr S z - 2H rz s) - ^e 2r ' '" s H z 2 E^ 

+e 2r2s r 2 (-z 4 E z2 H 2 z + z 2 B^H rr - 2z 2 E r E z H rr H rz - E^ H 2 r ) - r 2 A r H-}~s r 
+2r 2 H- 1 S r Z r + r 2 Y 2 - 2A r H- 1 S + 2A r H- 1 Z r + AH^SZr - 20Y 
+H- 1 (-r 6 z 2 H- 1 H 2 z S 2 r - Ar 4 z 2 H-^H 2 rz s~s r + 2r 4 z 2 H' 1 H 2 z ~s r Z r 

+2r 4 z 2 H rz s r s z - 2r 2 z 2 A r H- r 1 H 2 z s - 2r 2 z 2 X r Z H rz Y - r 2 z 2 H rr s 2 z 
-Ar 2 z 2 H- r l H 2 . z s 2 + Ar 2 z 2 H- 1 Hl z sZ r - 3r 2 z 2 H~ r 1 H 2 z s r 
+Ar 2 z 2 H rz ss z — 2r 2 z 2 H rz s r Z z — 2r 2 z 2 H rz s z Z r — r 2 H rz s r 

-.2 2 fj - , o~2-2 , o„2£V - fy c^2£V-l£V2 - 



+2z z A z H rz s + 2z z x z rz + 2z z H rr s z Z z - 6z z H- L H z J 
-Az 2 H rz sZ z + Az 2 H rz s z + H rr s z + 2H rz s) 
+H- 2 z 2 H 2 z (r 4 z 2 H- 1 H 2 J r - r 4 H rz s r - 2r 2 z 2 H i : r 1 H 2 z S - r 2 z 2 H rz ~s z 
+r 2 H rr S z + 2r 2 H rz l - z 2 H-^H 2 rz ) . 

Note that the expressions are manifestly regular on the axes and even in r and z. 
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APPENDIX D: MATTER EVOLUTION EQUATIONS 
1. Conservation form 



The matter evolution equations (|16H19[) can clearly be written in conservation form 

d t u+ [aF D (u)] tD = aS(u). 
Following , we replace pu with p^ = pu — a (kinetic energy) and regard as the set of conserved variables 

u= ( PK , J A , 3*, of. (Dl) 

The fluxes are 



and the source terms arc 



•> Pi< 


= J D -H D 


J~ Ja 


— Sa d , 


J- Jf 


= s D , 


•J a 


= £ D , 







S PK = (£ A -J A )(L> 7 - i ■ r - ' ' ■ * AJJ r i 

+\ 2 E A S A , 

Sj A = S AB (A B + L B ) + J A ( X + K v *) - A APH + L A r 

+\ 2 (E A J* + e AB S B B^), 
S JV = -(D I A + 3L A )S A + J^( X + SK^), 
S a = -(D 1 A + L A )V A + a( X + K v *) . 

2. Perfect fluid 

To evaluate the characteristic structure, we need to specify the matter model. Here, we consider a perfect fluid 
with four- velocity u a , normalized such that 

u a u a = —1 , 

rest mass density p, pressure p and internal energy e. The dependence of the pressure on the density and the internal 
energy is given by the equation of state 

p = p(p,e). (D2) 

With those definitions, the number density is 

N a = pu a 

and the energy-momentum tensor is given by 

T af3 = phu a vP + pg af3 , 

where h is the specific enthalpy, 



P 

The Lorentz factor is defined as 



h = l + e + -. (D3) 



W = -u a n a . (D4) 
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Observers who are at rest in a slice £(t) (i.e., who have four-velocity n a ) measure a coordinate velocity 

v A = W- l h a A u a , 



and the angular velocity is 



Hence we obtain the familiar relation 



where 



The variables 



W = {l-v 2 )- 1 ' 2 , (D5) 
v 2 = v A v A + X 2 v ip2 . 



w = (vA,vf,p,p,e, h,W) T 

are called primitive variables. 6 The conserved variables can be expressed in terms of the primitive variables as 

p K = phW 2 — p — pW , 
J a = phW 2 VA , 

Jf _ p hW 2 v ip , (D6) 
a = pW , 

and the remaining matter variables are 

t = phW 2 X 2 v ip2 + p, 
S A = phW 2 v v v Al 
S AB = phW 2 v A v B + P H AB , (D7) 
Y>A = pWv A . 

3. Characteristic decomposition 

The characteristic decomposition for 3+1 general relativistic hydrodynamics was first derived by the Valencia 
group qJ. The application to our (2+l)+l system is straightforward. Our method differs slightly in that we choose a 
general orthonormal basis (p A j J[ A ) in two-space as in section fWI and project vectors along p (index ||) and n (index 
_L). Following the notation of Q, we introduce a few abbreviations. From the equation of state l)D2(l . we form 

dp dp 2 _ p 

op oe p z 

where c s is known as the sound speed. Also set 7 



1 — Vi\Xt ' 1 — VnXf ' 

«I|-V ± , f=l-«jf, /\^h 3 W{l-JC- 1 ){C+ -C~)Z. 



6 Note only five of those are independent because of 1D2L 1D31 and 1D5I . 

7 Our definitions of £ and A differ from those in by a factor of A 2 to ensure regularity (see section lVllll . We have defined /C _1 instead 
of K, to allow for the special case of the ultrarclativistic equation of state, for which K~ 1 = 0. As a consequence, A above has been 
multiplied by K~ 1 and the characteristic variable £q,i has been divided by K.~ 1 . 
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The system is found to be strongly hyperbolic. The characteristic speeds in the /^-direction are 

Ao = v\\ , 

1 



A± = 



1 — v 2 cj 



v [{ (1 - <? s ) ± c s J (1 - «a) (1 - v*cl) - «f (1 - <*) 



The characteristic variables (corresponding to the left eigenvectors) are 



^0,1 
^0,2 
^0,3 



— j — U_l(<T + PA') + V||U_L J|| + (1 - ) Jj_| , 



^{-^(tr + pjcJ+^wyJy +(1-^)^} , 



ft 2 

— {r'/iP^u + [/cr 1 - ^ - (2 - jc -1 )^] J,, 

+ (2 - /C- 1 )V ± W r2 C(w|| J|| + v± J± + X 2 v v J ip ) 

+ [(JC- 1 - 1) (-«|| + V ± {W 2 i - 1)) - W 2 V±(\ (a + p K )} 

The inverse transformation (corresponding to the right eigenvectors) is given by 



J\\ 
J± 
J* 

PK 



-^-Mi + W(«iZo,2 + A 2 w%, 3 ) + 1+ + 17 , 
nW 

K.-\\l 0t i + 2hW 2 v ll (v ± l , 2 + A 2 v%, 3 ) + hW(C + lf + C~l~) , 
K~Vio,i + Wo,a + 2/iWV(uJo,2 + AVJo,3) + hWv±(l+ + l~) , 
K- l v v lo,i + W ,3 + 2WV(v_Li ,2 + AVl ,3) + hWv v (lf + 

+^M+i+ + ^-z.-)-z+-i7. 



4. Transformation from conserved to primitive variables 

The conserved matter variables (|D 1|) are the ones that are evolved in a numerical algorithm. To compute the 
remaining matter variables l|D7|l and the eigenvectors, the primitive variables have to be calculated from the conserved 
variables as an intermediate step. This transformation is much more involved than the opposite direction (|D6|) . To 
make it explicit, we have to specify an equation of state. Here, we consider the ultrarelativistic equation of state, 

p = (r - l)p tot = (r - l)p(e + 1) = ^r^ph , (D8) 

where p tot is the total energy density. 

Suppose we are given the conserved variables, and also form pjj — pk + cr. Consider the quantity 

J 2 = J A J A + X 2 J v2 . 

Using l|D6|) . I|D8|) and JD5J, we can express J 2 and pn in terms of the primitive variables as 

2 

72 I 1 \ „2ti/2/ttt2 



J = ^f^rj P 2 w 2 (w 2 -i) 

Ph = P (y&jW 2 - 1 
Eliminating W yields an equation for the pressure in terms of conserved variables: 



p = -2f3 PH + jA(3 2 p 2 H + (T - l)( P 2 H - J 2 ) 
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where /3 = (2 - T) /4. Next define 



_ (r-i)j A v _ (r - 1) 2 _ , a v2 

lp Lp 



We identify x A = W 2 v A and x v — W 2 v v and hence with l|D5|l we obtain 

f 

This now enables us to calculate the velocities, 

va = W~ 2 xa, v^ = W- 2 X v . 

The form of W~ 2 in l|D9l) guarantees that \va\, \v v \ < 1. This is most important since evolved speeds greater than 
unity (i.e. greater than the speed of light) can easily cause the numerical code to crash [Tfil |. 

Finally, we can calculate the specific enthalpy and rest mass energy density from (|D6|) and (|D8|> . 

h = — ^777 > P = 



avfW ' (r-i)/i 
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